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Abstract — Real number calculations on elementary functions 
are remarkably difficult to handle in mechanical proofs. In this 
paper, we show how these calculations can be performed within 
a theorem prover or proof assistant in a convenient and highly 
automated as well as interactive way. First, we formally establish 
upper and lower bounds for elementary functions. Then, based 
on these bounds, we develop a rational interval arithmetic where 
real number calculations take place in an algebraic setting. In 
order to reduce the dependency effect of interval arithmetic, 
we integrate two techniques: interval splitting and taylor series 
expansions. This pragmatic approach has been developed, and 
formally verified, in a theorem prover. The formal development 
also includes a set of customizable strategies to automate proofs 
involving explicit calculations over real numbers. Our ultimate 
goal is to provide guaranteed proofs of numerical properties with 
minimal human theorem-prover interaction. 

Index Terms — Real number calculations, interval arithmetic, 
proof checking, theorem proving 



I. Introduction 

Deadly and disastrous failures [l]-[3] confirm the shared 
belief that traditional testing, simulation, and peer-review are 
not sufficient to guarantee the correctness of critical soft- 
ware. Formal Methods in computer science refers to a set 
of mathematical techniques and tools to verify safety prop- 
erties of a system design and its implementation functional 
requirements. In the verification of engineering applications, 
such as aerospace systems, it is often necessary to perform 
explicit calculations with non-algebraic functions. Despite 
all of the developments concerning real analysis in theorem 
provers [4]-[8], the formal verification of the correctness of 
these calculations is not routine. 

Take, for example, the formula 
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180 v v 180 ' ~ 180 
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where g is the gravitational force and v = 250 knots is 
the ground speed of an aircraft. This formula appears in 
the verification of NASA's Airborne Information for Lateral 
Spacing (AILS) algorithm [9]. It states that the turn rate of 
an aircraft flying at ground speed v with a bank angle of 
35° is about 3° per second. A direct proof of this formula is 
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about a page long and requires the use of several trigonometric 
properties. 

In many cases the formal checking of numerical calculations 
is so cumbersome that the effort seems futile; it is then 
tempting to perform the calculations out of the system, and 
introduce the results as axioms. 1 However, chances are that 
the external calculations will be performed using floating-point 
arithmetic. Without formal checking of the results, we will 
never be sure of the correctness of the calculations. 

In this paper we present a set of interactive tools to automat- 
ically prove numerical properties, such as Formula (1), within 
a proof assistant. The point of departure is a collection of lower 
and upper bounds for rational and non-rational operations. 
Based on provable properties of these bounds, we develop a 
rational interval arithmetic which is amenable to automation. 
The series approximations and interval arithmetic presented 
here are well-known. However, to our knowledge, this is the 
most complete formalization in a theorem prover of interval 
arithmetic that includes non-algebraic functions. 

Our ultimate goal is to provide guaranteed formal proofs 
of numerical properties with minimum human effort. As 
automated processes are bound to fail on degenerate cases and 
waste time and memory on simple ones, we have designed a 
set of highly customizable proof strategies. The default values 
of the parameters are sufficient in most simple cases. However, 
a domain expert can set these parameters to obtain a desired 
result, e.g., the accuracy of a particular calculation. 

This paper merges and extends the results presented in 
[10], [11]. The rest of this document is organized as follows. 
Section II defines bounds for elementary functions. Section III 
presents a rational interval arithmetic based on these bounds. 
Section IV describes a method to prove numerical proposi- 
tions. The implementation of this method in a theorem prover 
is described in Section V. Last section summarizes our work 
and compares it to related work. 

The mathematical development presented in this paper has 
been written and fully verified in the Prototype Verification 
System (PVS) [12] 2 . PVS provides a strongly typed specifica- 
tion language and a theorem prover for higher-order logic. 
It is developed by SRI International. Our development is 
freely available on the Internet. The results on upper and 
lower bounds have been integrated to the NASA Langley PVS 
Libraries 3 and the rational interval arithmetic and the PVS 
strategies for numerical propositions are available from one of 

'As a matter of fact, the original verification of NASA's AILS algorithm 
contained several such axioms. 

2 PVS is available from http://pvs.csl.sri.com. 

3 http : / / shemesh .larc.nasa. gov/ f m/ ftp/larc/ 
PVS-library/pvslib . html. 



the authors 4 . 

For readability, we will use standard mathematical notations 
along this paper and PVS notations will be limited to illustrate 
the use of the library. In the following, we use the first letters 
of the alphabet a,b, . . . to denote rational numbers, and the 
last letters of the alphabet . . .x,y, z to denote arbitrary real 
variables. We use boldface for interval variables. Furthermore, 
if x is an interval variable, x denotes its lower bound and x 
denotes its upper bound. 



II. Bounds for Elementary Functions 

A PVS basic theory of bounds for square root and trigono- 
metric functions was originally proposed for the verification 
of NASA's AILS algorithm [9]. We have completed it and 
extended with bounds for natural logarithm, exponential, and 
arctangent. The basic idea is to provide for each real function 
/ : 1 h R, functions / : (R, N) h-> R and / : (R, N) h-> R 
closed under Q, such that for all x, n 



/(x,n) < f(x) 

l(x,n) < 

7(x,n + l) < 

lim f(x,n) = f(x) 



< f(x,n), 
f(x,n + l), 

= lim f{x,n) 



(2) 
(3) 
(4) 
(5) 



Formula (2) states that / and / are, respectively, lower and 
upper bounds of /, and formulas (3), (4), and (5) state that 
these bounds can ultimately be improved, as much as needed, 
by increasing the approximation parameter n. 

For transcendental functions, we use taylor approximation 
series. We performed a coarse range reduction [13] since the 
convergence of taylor series is usually best for small values. 
More elaborate range reduction techniques [14] would signif- 
icantly enhance the speed and the accuracy of the functions 
defined in Sections II and III. All the stated propositions in this 
section have been formally verified in the verification system 
PVS. 



A. Square root 

For square root, we use a simple approximation by Newton's 
method. For x > 0, 

sqrt(x, 0) = x + 1, 

1 x 

sqrtfx, n + 1) = — (y H — ), where y = sqrtfx, n), 

2 y 

sqrtfx, n) = == -. 

sqrtfx, n) 

Proposition 1: Vx > 0,n : < sqrt(x, n) < \fx < 
sqrt(x, n). 

The first inequality is strict when x > 0. 

4 http: / /research. nianet . or g/~munoz /Interval. 



B. Trigonometric functions 

We use the partial approximation by series. 

m 2i— 1 

sin(x,Ti) = E^ 1 )" 1 ^!) 
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m+l 2i— 1 

Mx,n) = E(- 1 )" 1 ( 2 ^_ 1)r 

m+1 x 2i 
cos(x,n) = 1 + ^(-1)' ' 



(2i)V 



m x 2i 
cos(x,n) = 1 + XI( _1 )'(27)I 



where m = 2n if x < 0; otherwise, m = 2n + 1, 

Proposition 2: Vx,n: sin(x,n) < sin(x) < sin(x,n). 
Proposition 3: Vx, n: cos(x,n) < cos(x) < cos(x, n). 

C. Arctangent and ir 

We first use the alternating partial approximation by series 
for < x < 1. 

atan fx, n) = > ' x Zl+L ^, ~\ , if < x < 1, 



2i+l' 



2n+l . 

i=i 

2n (—lY 

atanfx, n) = Vi 2l+1 - — 

y 1 ' ^ 2i + 1 

4=1 



, if0<x<l. 



We note that for x = 1 (which we might naively wish to use to 

define ir/4 and hence ir) the series: 1 — 5 + 5 — \~ >! ~\ ' does 

converge, but very slowly. Instead, we use the equality j = 
4 atan(l/5) — atan(l/239), that has much better convergence 
properties. Using this identity we can define bounds on ir: 



7r(n) = 16 atan fl, n) — 4 atanfl, n), 



7r(n) = 16 atanfl, n) — 4 atan fl, n). 

Proposition 4: Vn : Tr(n) < it < W(n). 
Now, using properties of arctangent, we extend the range of 
the function to the whole set of real numbers: 



atan fO, n) = atan(0, n) = 0, 



atan fx, n) 



7r(n) 



atanf — , n), if 1 < x, 
2 x 



atan fx, n) = — atan(— x,n), if x < 0, 

7f(n) 1 

atanfx, n) = — atan f — , n), if 1 < x, 

2 x 



atanfx, ri) 



- atan f— x, n), if x < 0. 



Proposition 5: Vx, n : atan fx, n) < atan(x) < 



atan(x, n). 

These are strict inequalities except when x = 0. 

The PVS definition of bounds on atan and 7r are presented 
in Listing 1. PVS developments are organized in theories, 
which are collections of mathematical and logical objects 
such as function definitions, variable declarations, axioms, 
and lemmas. The atan_approx theory first imports the 



definition of the arctangent function. Then, it declares variables 
n,x,px of types nat (natural numbers), real (real num- 
bers), and posreal (positive real numbers), respectively. For 
the scope of the theory, these variables are implicitly univer- 
sally quantified. Though writing definitions, lemmas, theorems 
and specially proofs in PVS requires some training, reading 
theories is possible to anybody with a minimal background in 
logic. 



PVS Listing 1 Definition of bounds on atan and n 

atan_approx: THEORY 
BEGIN 



D. Exponential 

The series we use for the exponential function is 



We could directly find bounds for negative x from this series 
as, in this case, the series is alternating. However, we will 
subsequently find that it is convenient to show that our bounds 
for the exponential function are strictly positive, and this is not 
true for all x < 0. Yet, this property holds for — 1 < x < 0. 
We define 

2(n+l)+l i 

exp(x,n) = > — , if — 1 < x < 0, 

i=0 
2(n+l) i 

cxp(x,n) = ^ — , if — 1 < x < 0. 
i=a l ' 

Using properties of the exponential function, we obtain 
bounds for the whole set of real numbers: 

cxp(0,n) = cxp(0, n) = 1, 

, , x -W 
exp(x,n) = expf — : — -,n) , if x < — 1 


exp(x,n) = -, if x > 0, 

exp(— x,n) 

_ _ x -N 

exp(x,n) = exp( — : — -,n) , if x < — 1 
-\x\ 



exp(a;, n) 



1 



cxp(— x, n) 



, if x > 0. 



Notice that unless we can ensure that all of the bounding 
functions are strictly positive we will run into type-checking 
problems using the bound definitions for x > 0, e.g., 
l/cxp(— x, n) is only defined provided exp(— x,n) ^ 0. 

Proposition 6: Vx, n : < exp(x, n) < exp(a;) < 
cxp(ir, n). 

These are strict inequalities except when x = 0. 



IMPORTING atan 

n : VAR nat 
x: VAR real 
px : VAR posreal 

atan_pos_lel_ub (n, x) : real = 
atan_ser ies_n (x, 2+n) 

atan_pos_lel_lb (n, x) : real = 
atan_ser ies_n (x, 2+n+l) 

atan_pos_lel_bounds : LEMMA 
< x AND x <= 1 IMPLIES 

atan_pos_lel_lb (n, x) < atan(x) AND 
atan (x) < atan_pos„lel_ub (n, x) 

pi_lbn (n) : posreal = 

4* (4*atan_pos„lel„lb (n, 1/5) - 
atan_pos_lel„ub (n, 1/239) ) 

pi„ubn(n) : posreal = 

4* (4*atan_pos_lel„ub (n, 1/5) - 
atan_pos_lel_lb (n, 1/239) ) 

pi„bounds : THEOREM 

pi_lbn (n) < pi AND pi < pi_ubn (n) 

atan_pos_lb (n, px) : real = 
IF px <= 1 THEN 

atan_pos_lel_lb (n, px) 
ELSE 

pi_lbn (n) 12 - atan_pos_lel_ub (n, 1/px) 
ENDIF 

atan_pos_ub (n, px) : real = 
IF px <= 1 THEN 

atan_pos_lel_ub (n, px) 
ELSE 

pi_ubn(n)/2 - atan_pos_lel_lb (n, 1/px) 
ENDIF 

atan_lb(x,n) : real = 

IF x > THEN atan_pos_lb (n, x) 

ELSIF x = THEN 

ELSE -atan_pos_ub (n, -x) ENDIF 

atan_ub(x,n) : real = 

IF x > THEN atan_pos_ub (n, x) 

ELSIF x = THEN 

ELSE -atan_pos_lb (n, -x) ENDIF 

atan_bounds: THEOREM 

atan_lb(x,n) <= atan (x) AND 
atan (x) <= atan_ub(x,n) 



E. Natural Logarithm 

For < x < 1, we use the alternating series for natural 
logarithm: 



END atan_approx 



In(x + 1) = 



Therefore, we define 

2» 



hi(x,n) = ^ _ ; ifl<a;<2, 



i=l 
2n+l 



Hx,n) = £(-1) 



-i (x - iy 

i 

-i (*-!)' 



, if 1 < x < 2. 



i=i 



Using properties of the natural logarithm function, we obtain 

ln(l,n) = In(l,n) = 

]n(x,n) = — ln( — ,n), if < x < 1, 
x 

ln(x,n) = — ln( — ,n), if < X < 1. 



Finally, we extend the range to the whole set of positive reals. 
If x > 2, we find a natural number m and real number y such 
that x = 2 m y and 1 < y < 2, by using the following recursive 
algorithm similar in spirit to Euclidean division: 

lnnat (x : posreal , k : posnat ) : [ nat , posreal ] = 
if x < k then (0,x) 
else 

let (m,y) = lnnat (x/k,k) in 
(m+l,y) 
endif 

We next prove the following property: 

Proposition 7: \/x > 1, k > 1 : k m < x < k m+1 ,y < 
k,x = k m y, where (m,y) = lnnat (x, k). 

If (m,y) = lnnat(2,x), we observe that 



ln(a;) = ln(2 m y) = m ln(2) + m(y). 



Hence, 



ln(a;,n) = m ln(2, n) + ln(y, n), if x > 2, 
ln(a;,n) = m ln(2, n) + ln(y, n), if x > 2. 

Proposition 8: \/x>0,n: ln(x,n) < ln(a;) < ln(x,n). 
These are strict inequalities except when x = 1. 

III. Rational Interval Arithmetic 

Interval arithmetic has been used for decades as a standard 
tool for numerical analysis on engineering applications [15], 
[16]. In interval arithmetic, operations are evaluated on range 
of numbers rather than on real numbers. A (closed) interval 
[a, b] is the set of real numbers between a and 6, i.e., 



[a,b] 



{x | a < x < b}. 



The bounds a and b are called the lower bound and upper 
bound of [a, b], respectively. Note that if a > b, the interval 
is the empty set. The notation [a] abbreviates the point-wise 
interval [a, a]. 

Interval computations can be performed on the endpoints 
or on the center and the radius. For this work, we decided to 
work on rational endpoints. Trigonometric and transcendental 
functions for interval arithmetic are defined using the bounds 
presented in Section II. 

Listing 2 shows a few definitions from the PVS the- 
ory Interval. Dots are used to simplify the presentation 



and hide some technical parts. The theory defines the type 
Interval as a record with fields ub and lb of type rat 
(rational numbers), variables x,y of type real, variable n 
of type nat, and variables X, Y of type Interval. 

PVS Listing 2 Definition of interval arithmetic 



Interval 
BEGIN 

Interval 



THEORY 



x, y 
n 

X, Y 



TYPE = [# 

lb : rat, 

ub : rat 

#] 



VAR real 

VAR nat 

VAR Interval 



+ (X, Y) : Interval 
- (X, Y) : Interval 
-(X) : Interval 



[ | lb (X) +lb (Y) , 

ub (X) +ub (Y) | ] 
[ |lb(X)-ub(Y) , 

ub (X) -lb (Y) | ] 
[ l-ub(X) , 
-lb(X) ) | ] 
* (X, Y) : Interval = . . . 
/(X,Y): Interval = X * [ | 1/ub (Y) , 

l/lb(Y) | ] 



Abs (X) 
Sq(X) 
A (X,n) 

U (X, Y) 



Interval = 
Interval = 
Interval = 

Interval 



[ | min (lb (X) , lb (Y) ) , 
max (ub (X) , ub (Y) ) | 



END Interval 



If X is a PVS interval, lb (X) is the lower bound and 
ub (X) is the upper bound of X. In PVS, we define the 
syntactic sugar [|x,y|] to represent the interval [x,y]. 
Interval union x U y, written in PVS X U Y, is defined as 
the smallest rational interval that contains both x and y. 

The four basic interval operations are defined as fol- 
lows [17]: 

x + y = [x + y,x + y], 
x-y = [x-y,x-y], 

x x y = [min{xy, xy, xy, xy}, max{xy, xy, xy, xy}] , 
x/y = xx [=,-], ifyy>0. 

y y - 

We also define the unary negation, absolute value, and power 
operators for intervals: 

-x = [-x, -x], 

|x| = [min{|x|, |x|}, max{|x|, |x|}], if xx > 0. 
|x| = [0,max{|x|, |x|}], if xx < 0. 

[1] if n = 0, 

[x",x n ] if x > or odd?(n), 

[x n ,x n ] ifx<0 andeven?(n), 

[0, maxjx™, x™}] otherwise. 



Interval operations are defined such that they include the 
result of their corresponding real operations. This property is 
called the inclusion property. 

Proposition 9 (Inclusion Property for Basic Operators): If 
x e x and y £ y then x®y £ x<X>y, where ® £ {+, — , x, /}. 
Moreover, — x £ — x, £ |x|, and x n £ x n , for n > 0. It 
is assumed that y does not contain in the case of interval 
division. 

Listing 3 specifies this property in PVS. The proposition x £ x 
is written x ## X. 

PVS Listing 3 Basic inclusion properties 

Add_inclusion : LEMMA 

x ## X AND y ## Y => x+y ## X+Y 

Sub_inclusion : LEMMA 

x ## X AND y ## Y =>• x-y ## X-Y 

Neg_inclusion : LEMMA 
x ## X -x ## -X 

Mult_inclusion : LEMMA 

x ## X AND y ## Y => x*y ## X*Y 

Div_inclusion : LEMMA 
NOT ## Y AND 

x ## X AND y ## Y ^> x/y ## X/Y 

Abs_inclusion : LEMMA 

x ## X => abs(x) ## abs(X) 

Sq_inclusion : LEMMA 

x ## X => sq(x) ## sq(X) 

Pow_inclusion : LEMMA 
x ## X => x A n ## X A n 



The inclusion property is fundamental to interval arithmetic. 
It guarantees that evaluations of an expression using interval 
arithmetic bound its exact real value. Any operation in interval 
arithmetic must satisfy the inclusion property with respect to 
its corresponding real operation. 

A. Interval comparisons 

There are several possible ways to compare intervals [18]. 
In this work, we use interval-rational comparisons and interval 
inclusions. 

x < a if x < a, similarly for <, 
x > a if x > a, similarly for >, 

x Q y if y < 2£ ar, d x < y. 

Proposition 10: Assume that x £ x, 

1) if x txi a then x x a, for X £ {<, <, >, >}, and 

2) if x C y then x £ y. 

We use ifa to denote >, >, <, or <, when XI is, respectively, 
<, <, >, or >. 

Proposition 11: If x XI a and x a, then x is empty. 
Notice that -i(x X a) does not imply x a. For instance, 
[—1,1] is neither greater nor less than 0. 



B. Square root, arctangent, exponential, and natural loga- 
rithm 

Interval functions for square root, arctangent, tt, exponen- 
tial, and natural logarithm are defined for an approximation 
parameter n > 0: 



[atan(x)] n 

Mr. 

[exp(x)]„ 
Pn(x)]„ 



[sqrt(x, n),sqrt(x, n)], if x > 0, 



[atan(x, n), atan(x, n)], 

[7r(n),7f(n)], 

[cxp(x, n), cxp(x, n)], 

[ln(x, n) , ln(x, n)], if x > 0. 



As consequence of Propositions 1,5, 6, and 8 in Section II, 
and the fact that these functions are increasing, the above 
functions satisfy the following inclusion property. 

Proposition 12: For all n, if x £ x then f(x) £ [/(x)]„, 
where / £ {«/", at an, exp, In}. Moreover, it £ [n] n . It is 
assumed that x is non-negative in the case of square root, 
and x is positive in the case of natural logarithm. 

C. Trigonometric functions 

Parametric functions for interval trigonometric functions are 
defined by cases analysis on quadrants where the functions 
are increasing or decreasing. The mathematical definitions are 
presented in Figure 1. 

Note that sin and cos are defined for the whole real line. 
However, for angles a such that \a\ > tt_ both functions will 
return the interval [—1,1], a valid bound but not a very good 
one. Furthermore, the expression n + 5 in Formula (8) is 
necessary to guarantee that lower and upper bounds of cosine 



are strictly positive in the interval 



7r(n+5) 
2 



7r(n+5) 



, and thus, 



the interval tangent function is always defined in that interval. 

The interval trigonometric functions satisfy the inclusion 
property. 

Proposition 13: If x £ x then f(x) £ [/(x)]„, where / £ 
{sin, cos}. Moreover, if x C [-2^1, a&^l], t aa(i) £ 
[tan(x)]„. 



The next section proposes a method to prove numerical 
propositions based on the interval arithmetic described here. 

IV. Mechanical Proofs of Numerical Propositions 

Arithmetic expressions are defined by the following gram- 
mar, where V is an denumerable set of real variables: 



a \ x \ e + e \ e — e | 
e/e | |e| | e 4 | Ve | 7T 
cos(e) | tan(e) | exp(e) 



-e | e x e | 
sin(e) | 
ln(e) | atan(e) 



Numerical propositions P have either the form e± \x e%, 
where x £ {<, <, >, >}, or the form e £ a, where a is a 
constant interval (an interval with constant rational endpoints). 
As usual, parentheses are used to group real and interval 
expressions as needed. 



a 


£ 


Q 


i 


£ 


N 


X 


e 


V 



[sin(x)]„ = 



[cos(x)], 



[tan(x)], 



Fig. 1. Interval trigonometric functions 



sin(x, n), sin(x, n)} 
sin(x, n), sin(x, n)} 
min{ sin(x, n) , sin(x, n) } , 1] 
-[sin(-x)]„ 
-1,1] 

cos(x, n), cos(x, n)} 
cos(-x)]„ 

min{cos(x, n), cos(x, n)}, 1] 



(") iWi 

2 ' 2 J 



xC[f,i(n) 



. sin , sin ,_ 

(x, n + 5 ) , (x, n 

cos COS 



if x C 
else if 
else if x C [0, 7r(n)], 
else if x C [— 7r(n),0] 
otherwise, 

if xC[0,7r(n)], 

else if x C [— 7r(n),0], 
else if 
otherwise 

7r(n + 5) 7r(n + 5) 



xc [_£fcl, sM] 



+ 5)], ifxC[- 



(6) 



(7) 



(8) 



A context T is a set of hypotheses of the form x G x. 
A ground context is a context where all the intervals are 
constant. In the following, we use logical judgments in the 
sequent calculus style, e.g., T h P, where all free variables 
occurring in P are in V. The intended semantics of a judgment 
r h Pis that the numerical proposition P is true under the 
hypotheses V. 

Given a context T, an approximation parameter n, and an 
expression e, such that the free variables of e are in T, we 
define the interval expression [e]£ by recursion on e. 



x, where (x G x) G T, 



[ei 



e2 



ei 



e 2 



where 



{+,-*,/}, 



= ([e} r n y, 



IW 



[/(Nn)]n, 

where / G {sin, cos, tan, exp, In, atan}. 



Theorem 1 (Inclusion): Let T be a context, n an approxi- 
mation parameter, and e a well-defined arithmetic expression 
in r, i.e., side conditions for division, square root, logarithm, 
and tangent are satisfied, 

r h e G [C (9) 
Proof: By structural induction on e and proposi- 
tions 4, 9, 12, and 13. ■ 

A. A general method for numerical propositions 

We propose a general method to prove numerical proposi- 
tions. First, consider a judgment of the form 

F h ei X e 2 , 

where T is a ground context. 

1) Select an approximation parameter n. 

2) Define e = e\ — e 2 . 

3) Evaluate [e]£ cxi 0. If it evaluates to true, the following 
judgment holds 



4) 



5) 



In that case go to step 5. 

Evaluate [e]£ tfo 0. If this evaluates to true then fail. By 
Proposition 11, the judgment T h [e]^ M cannot hold. 
If [e]Jj ^1 evaluates to false, increase the approximation 
parameter and return to step 3. 
By Theorem 1, 



r h e G [e]l 



6) Proposition 10 yields 



7) By definition, 



8) Therefore, 



T h etoO. 



T h ei - 62 CO 0. 



T h ei ixi e 2 . 



r h 



M o. 



The method above can be easily adapted to judgments of 
the form V h e IX a. In this case, the interval expression 
[e]n — a ' s evaluated. If the expression evaluates to true, then 
the original judgment holds by Theorem 1 and Proposition 10. 
Otherwise, the method should fail. 

The general method is sound, i.e., all the steps can be 
effectively computed and each one is formally justified. In 
particular, the propositions [e]„ M 0, [e]Jj t/l 0, and [e]£ C a 
can be mechanically computed as they only involve rational 
arithmetic and constant numerical values. The method is not 
complete as it does not necessarily terminate. Even if e only 
involves the four basic operations and no variables, it may be 
that both [e]£ M and [e]£ ^ evaluate to false. 

The absence of a completeness result is a fundamental 
limitation on any general computable arithmetic. At a practical 
level, the problem arises because all we have available are 
a sequence of approximations to the real numbers x and y; 
provided x and y differ, with luck we will eventually have 
a pair of approximations whose intervals do not overlap, and 
hence we can return a result for x cxi y. However, if x and y 
are the same real number (note we might not necessarily get 
the same sequence of approximations for both x and y), we 
can never be sure whether further evaluation might result in 
us being able to distinguish the numbers. 



B. Dependency effect 

The dependency effect is a well-known behavior of interval 
arithmetic due to the fact that interval identity is lost in interval 
evaluations. This may have surprising results, for instance 
x — x is [0] only if x is point-wise. Moreover, as we have 
seen in Section III-A, both x > a and x < a may be 
false. Additionally, interval arithmetic is subdistributive, i.e., 
xx (y + z ) Q xxy + xxz. In the general case the inclusion 
is strict and some dependency effects appear as soon as a 
variable is used more than once in an expression. 

For the method presented in Section IV- A, it means that the 
arrangement of the expression e matters. For instance, assume 
that we want to prove x E [0, 1] h 2 x x > x. This 
is pretty obvious in arithmetic as x is a non-negative real. 
Using our method, we first consider the arithmetic expression 
e = 2 x x — x and then construct the interval expression 
[e]n = 2 x x — x, where x = [0, 1]. For any approximation 
parameter n, [e]Jj evaluates to [—1, 2] which is neither greater 
nor less than 0. Therefore, the method will not terminate. On 
the other hand, if instead of the arithmetic expression 2 x x— x, 
we consider the equivalent arithmetic expression x, we have 
[x]Jj = [0, 1] and [0, 1] > evaluates to true. 

A second observation is that because of the dependency 
effect the width of intervals also matters. Consider again the 
expression e = 2 x x — x. We have seen that the interval 
evaluation of [e]£, for x G [0,1], results in [—1,2], which 
is not sufficient to prove that [e]£ > 0. On the other hand, 
the expression [e]£ evaluates to [—1/2,1] when x E [0,1/2] 
and it evaluates to [0,3/2] when x E [1/2,1]. Therefore, we 
can prove that, for x E [0,1], [e] r n C [-1/2,1] U [0,3/2], 
i.e., [e]£ C [—1/2,3/2], which is a better approximation than 
[—1,2]. If we continue dividing the interval [0,1] and com- 
puting the union of the resulting intervals, we can eventually 
prove that [e]£ + e > for an arbitrary small e > 0. 

These observations lead to two enhancements of the general 
method. First, we could divide each interval in T before 
applying the general technique. Second, we may want to 
replace the original expression by an equivalent one that is 
less prone to the dependency effect. 

C. Interval splitting 

In interval arithmetic, the dependency effect of the union 
of the parts is less than the dependency effect of the whole. 
Indeed, the simplest way to reduce the dependency effect is 
to divide the interval variables into several tiles (subintervals) 
and to evaluate the original expression on these tiles separately. 
This technique is called interval splitting or sub-paving and is 
expressed by the following deduction rule. 

Proposition 14: Let T be a context, e an expression whose 
free variables are x and those in T, e an interval expression, 
and x, Xi, . . . , x„ intervals such that x = Ui<i< n x «' 



V 1 < i < 77 : i£Xi,r h e E e 



[Splitting] 



i£x,r h e E e 
The Splitting rule can be iterated to obtain a splitting for 
multiple variables. Note that the number of tiles generated by 
interval splitting is exponential in the number of variables. 
Indeed, if ki is the number of tiles of the first variable alone, 



k 2 is the number of tiles of the second variables alone, and so 
forth, the total number of tiles to be considered for m variables 

is Ul<j< m k j- 

The integration of the Splitting rule into the general method 
is straightforward. First, a splitting is computed for a given set 
of variables in T. Then, the general method is applied to all 
cases. If the general method is successful in all of them, by 
Proposition 14, the original judgment holds. Otherwise, the 
method fails and a new splitting may be considered. 

D. Taylor Series Expansions 

Replacing 2 x x — x by x can be done automatically. In fact, 
as we will see in Section V, these kinds of simplifications are 
performed by our PVS implementation of the general method. 
However, these simplifications may not be sufficient even for 
simple expressions such as x x (1 — x), where x E [0.1]. The 
subdistributivity property of interval arithmetic states that the 
interval evaluation of x x (1 — x ) is better than that of the 
equivalent expression x — x 2 . Unfortunately, that evaluation 
is not good enough to prove that x X (1 — x) E [0, 1/4]. In 
this case, as a domain expert knows, the optimal answer is 
obtained with the equivalent expression 1/4— (1/2 — x) 2 . The 
solution is a lot less intuitive when non-algebraic functions are 
involved. 

Taylor's theorem states that a n-differentiable function can 
be approximated near a given point by a polynomial of degree 
?? whose coefficients depend on the derivatives of the function 
at that point. In interval arithmetic, taylor's theorem can be 
expressed by the following deduction rule. 

Proposition 15: Let x, xo, . . . , x„ be strictly proper inter- 
vals, / a n-differentiable function on a variable x E x, and 
c E x a constant, 

V < i < n : h /«(c) 6 ^ 
xex h /(">(aQ gx„ 
x E x h f{x) E Er =0 (x, x (x - cY)/ t \ ay MJ 
The expression of Taylor's rule shows that interval x appears 
only once in each term of order i for i between 1 and 77 — 1 
preventing any dependency effect due to x in a term alone. 
The term of order 77 suffers some dependency effect as x also 
appears in the definition of x n . In most cases, n = 2 is used 
to cancel first order dependency effects as presented Listing 4. 
But in cases where the first derivatives nearly vanish or where 
the evaluation of the last derivative introduces significant 
dependency effects, we compute more terms to reach some 
better bounds. 

Using Taylor's rule require more work than the Splitting 
rule. In particular, we need to provide intervals Xo,...,x n 
and constant c that satisfy the hypotheses of the rule. For c we 
choose the middle point of x unless the user proposes another 
point. It follow immediately that c E x. For < i < n, we 
choose Xj = [/^(c)]„ and, by Theorem 1, we have /W(c) E 
x^ Finally, we choose x„ = [/(") (x)]^, where T is the context 
x E x. By Theorem 1, we have T h /W(x) <E x„. 

In order to prove the judgment x E x h /(x) E a, we 
consider the interval expression £™ =0 (xi x(x — c) l )/7!Ca for 
a given n. If it evaluates to true, then the original judgment 
holds by Taylor's rule and Proposition 10. If the evaluation 



returns false, the method fails and a higher expansion degree 
n may be considered. 

For better results, the evaluation of £™ =0 (xi x (x— c) l )/i\ C 
a can be performed using the splitting technique. Contrary to 
the approach described in [19], we do not have to generate a 
new taylor approximation for each tile. By using an interval- 
based taylor expansion, the same expression can be reused 
for all the tiles. One single global taylor expansion has to be 
validated, and the proofs for all the tiles simply consist in an 
interval evaluation of this expansion. We do not suffer from the 
taylor coefficients being irrational numbers, they are simply 
given by interval expressions involving rational functions. 
Relying on rational interval arithmetic leads to conceptually 
simpler proofs. 

Section V describes how the general method and its ex- 
tensions are implemented in the PVS theorem prover and 
illustrates the practical use of the library with a few examples. 

V. Verified Real Number Calculations in PVS 

The interval arithmetic presented in this paper has been 
developed as a PVS library called Interval. This library 
contains the specification of interval arithmetic described here 
and the formal proofs of its properties. We believe that a 
domain expert can use this library with a basic knowledge of 
theorem provers. Minimal PVS expertise is required as most of 
the technical burden of proving numerical properties is already 
implemented as proof strategies. 

A. Strategies 

The numerical strategy is the basic strategy that im- 
plements the general method and its extensions described in 
Section IV. For instance, Formula 1 can be specified in PVS 
as follows (comments in PVS start with the symbol % and 
extend to the end of the line): 

g : posreal = 9.8 %[m/s A 2] 
v : posreal = 250*0.514 %[m/s] 

tr35: LEMMA 

(g*tan (35*pi/180) /v) * 180/pi 
## [I 3, 3.1 I ] 

%|- tr35: PROOF (numerical) QED 

We emphasize that, in PVS, tan and pi are the real math- 
ematical function tan and constant it, respectively. Lemma 
tr35 is automatically discharged by the numerical strat- 
egy, which can be entered interactively or in batch mode, as 
in this case, via the ProofLite library developed by one of the 
authors [20]. 

Another example is the proof of the inequality 4.1.35 
in [13]: 

Vx: < x < 0.5828 =S> |ln(l-x)|<y. 
The key to prove this inequality is to prove that the function 

G(x) = y-ln(l-z) 



satisfies G(0.5828) > 0. In PVS: 

G(x|x < 1): real = 3*x/2 - ln(l-x) 

A_and_S : lemma G (0.5828) > 
%|- A_and_S : PROOF (numerical : clefs "G") QED 

In this case, the optional parameter :defs "G" tells 
numerical that the user-defined function G has to be 
expanded before performing the numerical evaluation. The 
original proof of this lemma in PVS required the manual 
expansion of 19 terms of the In series. 

The numerical strategy is aimed to practicality rather 
than completeness. In particular, it always terminate and it 
is configurable for better accuracy (at the expense of perfor- 
mance). 

Termination is trivially achieved as the strategy does not 
iterate for different approximations, i.e., step 3 either goes 
to step 5 or fails. In other words, if numerical does not 
succeed, it does nothing. Furthermore, numerical uses a 
default approximation parameter n = 3, which gives an 
accuracy of about 2 decimals for trigonometric functions. 
However, the user can increase this parameter or set a different 
approximation to each function according to his/her accuracy 
needs and availability of computational power. Currently, there 
is no direct relation between the approximation parameter 
and the accuracy, as all the bounding functions have different 
convergence rates. On-going work aims to provide, an absolute 
error of at most 2~ p for any expression with a new approx- 
imation parameter p. The strategy has not been designed to 
reuse past computations. Therefore, it will be prohibitively 
expensive to automatically iterate numerical to achieve a 
small approximation on a complex arithmetic expression. 

In order to reduce the dependency effect, the numerical 
strategy automatically rearranges arithmetic expressions using 
a simple factorization algorithm. Due to the subdistributivity 
property, the evaluation of factorized interval expressions is 
more accurate than that of non-factorized ones. A set of 
lemmas of the NASA Langley PVS Libraries are also used 
as rewriting rules on arithmetic expressions prior to numerical 
evaluations. This set of lemmas is parameterizable and can 
be extended by the user. For instance, trigonometric functions 
applied to notable angles are automatically rewritten to their 
exact value. Therefore, numerical is able to prove that 
sin(7r/2) 6 1, even if this proposition is not provable using our 
interval arithmetic operators alone. Although it is not currently 
implemented, this approach can also be used to normalize 
angles to the range [— it, ir] that is suitable for the interval 
trigonometric functions in Sections III-C. 

The splitting technique is implemented by allowing the 
user to specify the number of tiles to be considered for each 
interval variable or a default value for all of them. The strategy 
will evenly divide each interval. For example, the simple 
expression in Section IV-D can be proven to be in the range 
[0,9/32] using a splitting of 16 subintervals. 

fair : LEMMA 

x ## [ |0,1| ] IMPLIES x* (1-x) ## [ | 0,9/32 |] 



%|- fair : PROOF (instint : splitting 16) QED 

In this example we have used the instint strategy. 
This strategy is built on top of numerical and per- 
forms some basic logic manipulations such as introduction 
of real variables and interval constants. In this case, the 
proof command (initint : splitting 16) is equiv- 
alent to (then (skeep) (numerical :vars ("x" 
" [ I , 1 | ] " 1 6 ) ) ) . It instructs PVS to introduce the real 
variable x and then to apply numerical by splitting 16 times 
the interval [0, 1]. 

The taylor series expansion technique is implemented in 
two steps. First, the taylor strategy automatically proves 
Proposition 15 for a particular function / and degree n. In 
the following example, we show that x E x h x x (1 — x) 6 
Si=o( Xi x ( x — c )')/*!' provided that x is strictly proper. 

F(X) : MACRO Interval = X*(l-X) 
DF(X) : MACRO Interval = 1 - 2*X 
D2F(X): MACRO Interval = [| -2 |] 

ftaylor : LEMMA 

x ## X AND StrictlyProper? (X) IMPLIES 
x*(l-x) ## Taylor2 [X] (F, DF, D2F) 

%|- ftaylor : PROOF (taylor) QED 

The keyword MACRO tells the theorem prover to automati- 
cally expand the definition of the function. The expression 
Taylor2[X] (F , DF , D2F ) corresponds to X)Lo( Xj x ( x - 
c) l )/i\, where F, DF, and D2F are the interval functions 
corresponding to /, it 1st, and its 2nd derivative. 

Finally, the strategy instint is called with the lemma 
ftaylor. 

best : LEMMA 

x ## [ |0,1| ] IMPLIES x*(l-x) ## [|0, 1/4 | ] 

% | - best : PROOF 

%|- (instint :taylor "ftaylor") 

%|- QED 

B. A simple case study 

The arctangent function is heavily used in aeronautic appli- 
cations as it is fundamental to many Geodesic formulas 5 . One 
common implementation technique uses an approximation 
of the arctangent on the interval x = [—1/30,1/30] after 
argument reduction [21]. For efficiency reasons, one may want 
to approximate the function atan(x) to single precision by the 
polynomial 

. . 11184811 , 13421773 5 

fix) - x — x — X 

y ' " 33554432 67108864 

The coefficients of the polynomial approximation are stored 

exactly using IEEE single precision. 

The objective of this case study is to show that 

xe [-1/30,1/30] h atan(x)-r(x) G [-2^,2'% 

for different values of i. The PVS specification of this problem 
for some values of i is presented in Listing 4. All the lemmas 

5 See, for example, Ed William's Aviation Formulary at http:// 
williams -best . vwh . net /avf orm . htm. 
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Fig. 2. Time required to prove tan(a:) — r(x) £ [—1/30, 1/30] 

are automatically discharged by the instint strategy with 
different splitting and taylor's expansion degrees. As expected 
taylor's expansions and splitting get better results than splitting 
alone. Moreover, second degree expansions are almost always 
better than first degree expansions. This is not necessarily 
the case as illustrated by lemmas f air_atan_tl_14 and 
f air_atan_t2_14: for i = 14, a first degree expansion 
with no splitting is enough to prove the property, while a 
second degree expansion requires a splitting of 2. 

On a tile t of x, the width of the error expression E that 
does not use taylor's theorem evaluated on t is larger than 
the sum of the width of expressions At an and R. As the 
derivative of the arctangent is between 0.9989 and 1 on x, we 
could expect that the width of R is at least twice the width 
of tile t. Therefore, to obtain an error bound of [— 2~ l ,2~ l ] 
we cannot use tiles larger than 2~ J and we will need at least 
2 l /15 w 2 l - 1A tiles. 

We use the same kind of simple calculation to show that 
since \e'{x)\ < 2.37- 10~ 6 we will need about 2 l ~ 14 - 8 tiles of 
width 2~ l ■ 10 6 /2.37. This figures are accurate when we use 
second degree expansion but actual computations may require 
more tiles due to some dependency effects introduced when 
we use first degree expansions. 

Figure 2 presents a summary of the time required to prove 
tan(as) - r{x) G [-1/30,1/30] for i in the range [0,20] 
using splitting, splitting and first degree taylor's expansion, 
and splitting and second degree taylor's expansion. 

C. Implementation and Performance Issues 

Actual definitions in PVS have been slightly modified for 
efficiency reasons. For instance, multiplication is defined using 
a case analysis on the sign of the operands. Additionally, 
all interval operations are completed by returning an empty 
interval if side conditions are not satisfied. This technique 
avoids some type correctness checks that are expensive. 

The strategies in this library work over the PVS built-in 
real numbers. The major advantage of this approach is that the 
functionality of the strategies can be extended to handle user 
defined real functions without modifying the strategy code. 
Indeed, optional parameters to the numerical strategy allow 
for the specification of arbitrary real functions. If the interval 
interpretations are not provided, the strategy tries to build them 
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fair_atan : THEORY 
BEGIN 



x 

r (x) 
e (x) 
Xt 



var real 
MACRO real 
MACRO real 
Interval 



x - (11184811/33554432) 
atan (x) - r (x) 
[| -1/30, 1/30 |] 



x A 3 



(13421773/67108864) * x A 5 



fair_atan_8 : LEMMA x ## Xt IMPLIES e (x) ## [|-2 A -8, 2 A -8|] 
%|- fair_atan_8 : PROOF (instint : splitting 18) QED 



X : var Interval 

R(X) : MACRO Interval = X - 11184 811/33554432 * X A 3 - 13421773/67108864 * X A 5 

E(X) : MACRO Interval = Atan (X, 4) - R(X) 

DE(X) : MACRO Interval = 

1 / (1 + Sq(X)) - 1 + 3* (X A 2* (11184811/33554432) ) + 5 * (X A 4 * ( 1342 17 73/ 67 1 08 8 64 ) ) 



atan_taylorl : LEMMA StrictlyProper ? (X) AND x ## X IMPLIES e (x) ## Taylorl [X] (E, DE) 
%|- atan_taylorl : PROOF (taylor) QED 

f air_atan_tl_14 : LEMMA x ## Xt IMPLIES e (x) ## [|-2 A -14, 2 A -14 | ] 
%|- f air_atan_tl_14 : PROOF (instint : taylor " atan_taylor 1 " ) QED 

f air_atan_tl_2 : LEMMA x ## Xt IMPLIES e (x) ## [|-2 A -20, 2 A -20|] 
%|- f air_atan_tl_20 : PROOF (instint :taylor " atan_taylor 1 " : splitting 13) QED 

D2E(X) : MACRO Interval = 

-2*X/Sq(l + Sq(X)) + 2 * ( X A 3 * ( 1 34 2 1 7 7 3 / 67 1 8 8 6 4 ) ) + 6* ( ( 1 11 84 81 1/33554 432 ) *X) 



atan_taylor2 : LEMMA StrictlyProper? (X) AND x ## X IMPLIES e (x) ## Taylor2 [X] (E, DE, D2E) 
%|- atan_taylor2 : PROOF (taylor) QED 

f air_atan_t2_14 : LEMMA x ## Xt IMPLIES e (x) ## [|-2 A -14, 2 A -14 | ] 
%|- f air_atan_t2_14 : PROOF (instint :taylor " atan_taylor2 " :spitting 2) QED 

f air_atan_t2_2 : LEMMA x ## Xt IMPLIES e (x) ## [|-2 A -20, 2 A -20|] 
%|- f air_atan_t2_2 : PROOF (instint :taylor " atan_taylor2 " splitting 5) QED 



END fair_atan 



from the syntactic definition of the functions. The trade-off 
for the use of the PVS type real, in favor of a defined data 
type for arithmetic expressions, is that the function [e]^ and 
Theorem 1 are at the meta-level, i.e., they are not written in 
PVS. It also means that the soundness of our method cannot 
be proven in PVS itself. In particular, Theorem 1 has to be 
proven for each particular instance of e and [e]£. This is not 
a major drawback as, in addition to numerical, we have 
developed a strategy called inclusion that discharges the 
sequent T h e € [e]JJ whenever is needed. PVS strategies are 
conservative in the sense that they do not add inconsistencies 
to the theorem prover. Therefore, if numerical succeeds to 
discharge a particular goal the answer is correct. 

Finally, our method relies on explicit calculations to evaluate 
interval expressions. In theorem provers, explicit calculations 
usually means symbolic evaluations, which are extremely 
inefficient for the interval functions that we want to calculate. 
To avoid symbolic evaluations, numerical is implemented 
using computational reflection [22]-[24]. Interval expressions 
are translated to Common Lisp (the implementation language 
of PVS) and evaluated there. The extraction and evaluation 
mechanism is provided by the PVS ground evaluator [25]. 
The result of the evaluation is translated back to the PVS 
theorem prover using the PVSio library developed by one of 
the authors [26]. 



VI. Conclusion and Limits of Tractability 

We have presented a pragmatic approach to verify ordinary 
real number computations in theorem provers. To this end, 
bounds for non-algebraic functions were established based on 
provable properties of their approximation series. Furthermore, 
a library for interval arithmetic was developed. The library 
includes strategies that automatically discharges numerical 
inequalities and interval inclusions. 

The PVS Interval library contains 306 lemmas in total. It 
is roughly 10 thousand lines of specification and proofs and 1 
thousand lines of strategy definitions. These numbers do not 
take into account the bounding functions, which have been 
fully integrated to the NASA Langley PVS Libraries. It is 
difficult to estimate the human effort for this development 
as it has evolved over the years from an original axiomatic 
specification to a fully foundational set of theories. As far 
as we know, this is the most complete formalization within 
a theorem prover of an interval arithmetic that includes non- 
algebraic functions. 

Research on interval analysis and exact arithmetic is rich 
and abundant (see for example [17], [27], [28]). The goal 
of interval analysis is to compute an upper bound of the 
round-off error in a computation performed using floating- 
point numbers. In contrast, in an exact arithmetic framework, 
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Fig. 3. Alternate f air_atan theorems will make use of interval arithmetic 



an accuracy is specified at the beginning of the computation 
and the computation is performed in such way that the final 
result respects this accuracy. 

Real numbers and exact arithmetic is also a subject of 
increasing interest in the theorem proving community. Pio- 
neers in this area were Harrison and Gamboa who, indepen- 
dently, developed extensive formalizations of real numbers for 
HOL [4] and ACL2 [6]. In Coq, an axiomatic definition of 
reals is given in [7], and constructive definitions of reals are 
provided in [29] and [30]. As real numbers are built-in in 
PVS, there is not much meta-theoretical work on real num- 
bers. However, a PVS library of real analysis was originally 
developed by Dutertre [31] and currently being maintained 
and extended as part of the NASA Langley PVS Libraries. 
An alternative real analysis library is proposed in [8]. 

Closer to our approach are the tools presented in 
[32] and [10]. These tools generate bounds on the round-off 
errors of numerical programs, and formal proofs that these 
bounds are correct. The formal proofs are proof scripts that 
can be checked off-line using a proof assistant. 

Our approach is different from previous works in that we 
focus on automation and pragmatism. In simple words, our 
practical contribution is a correct pocket calculator for real 
number computations in formal proofs. Thanks to all the 
previous developments in theorem proving and real numbers, 
lemmas like Lemma tr35 and Lemma A_and_S are prov- 
able in HOL, ACL2, Coq, or PVS. The Interval library make 
these proofs routine in PVS. 

As in real life, users benefit in managing both a pocket 
calculator and a graphic tool. The fact that the example 
proposed in Section V-B is reaching the limits of tractability is 
not a problem. Our library aims at providing some simple tools 
that can be used seamlessly in proofs. Figure 3 would prompt a 
careful user that f air_atan theorems are a consequence of 
the fact that the derivative of the error is always positive. Such 
a fact could happen to be difficult to prove leading some one to 
prove that the error is bounded on some subintervals and that 
the derivative is always positive on some other subintervals. 
Anyways, such proofs will involve our library more than once. 

We continue developing this library and it is currently 
being used to check numerical properties of aircraft navigation 
algorithms developed at the National Institute of Aerospace 



(NIA) and NASA. Future enhancements include: 

• Development of a fully functional floating point arith- 
metic library [33] in order to generate guaranteed proofs 
of round-off-errors [32]. 

• Integration of this library and an exact arithmetic formal- 
ization in PVS developed by one of the authors [34]. 

• Implementation of latest developments on Taylor Mod- 
els [35]-[37], which will enable a greater automation of 
the Taylor's series expansion technique. 
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